# Author: Gregory Smith
# Date: 12 - 3 - 2018 
# Title: Smith 2018 Secret But Constrained Visualizations 
# Job: Includes the Code for Every Visualization Except the AME and Interaction Plots

#---- Packages 


list.of.packages <- c("tidyverse", "readxl", "scales", "cowplot", "kableExtra",
                       "janitor", "cshapes", "GISTools")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)


#--- Plots and Data Manipulation 
library(tidyverse)
library(readxl)
library(scales)
library(cowplot) # puts plots into a grid 
library(kableExtra) # Creates latex tables 
library(janitor)

#---- Maps 
library(cshapes)     
library(GISTools)    

# Change Working Directory to Replicate Code 
setwd("~/Dropbox/covert_ops_paper/2019-smith-io-secret-but-constrained-replication-files/")


#---- Data 
dta <- read_csv("data/smith_2019_secret_but_constrained_dta.csv")
us_domestic_data <- read_csv("data/final_us_domestic_year.csv")

bvp_onset_dta <- read_excel("smith-2019-bvp-marginal-effects/smith_2019_bvp_onset_marginal_effects.xlsx", trim_ws = TRUE)
bvp_ongoing_dta <- read_excel("smith-2019-bvp-marginal-effects/smith_2019_bvp_ongoing_marginal_effects.xlsx", trim_ws = TRUE)

#--------- Plots of Frequency of Covert Operations Per Year 

#---- Number of Distinct Operations Initiated Per Country
onset_count_by_country <- dta %>%
  dplyr::group_by(ccode_current, isocode) %>% 
  dplyr::summarise(orourke_onset = sum(orourke_onset, na.rm = TRUE),
                   pooled_onset = sum(us_covert_onset, na.rm = TRUE))

#---- Number of Operations Active and Initiated Per Year
year_counts <- dta %>%
  dplyr::group_by(year) %>% 
  dplyr::summarise(orourke_ops_active = sum(orourke_covert_intervention, na.rm = TRUE),
                   orourke_ops_initiated = sum(orourke_onset, na.rm = TRUE),
                   pooled_ops_active = sum(us_covert_operation, na.rm = TRUE),
                   pooled_ops_initiated = sum(us_covert_onset, na.rm = TRUE),
                   major_gml_initiated = sum(major_us_onset, na.rm = TRUE),
                   minor_gml_initiated = sum(minor_us_onset, na.rm = TRUE),
                   major_force_initiated = sum(major_force_onset, na.rm = TRUE),
                   minor_force_initiated = sum(minor_force_onset, na.rm = TRUE))


covert_frequency <- left_join(year_counts, us_domestic_data, by = "year") %>%
  dplyr::select(-X1) %>% 
  dplyr::mutate(
    president = case_when(
      truman == 1      ~ "truman",
      eisen  == 1      ~ "eisenhower",
      kennedy == 1     ~ "kennedy",
      johnson == 1     ~ "johnson",
      nixon == 1       ~ "nixon",
      ford == 1        ~ "ford",
      carter == 1      ~ "carter",
      reagan == 1      ~ "reagan",
      bush   == 1      ~ "bush",
      TRUE             ~ "other"
    )
  )


#----  Frequency Visualizations  

#----- Operations Active Per-Year for the O'Rourke and Pooled Samples 

orke_active <- ggplot(covert_frequency, aes(x = year, y = orourke_ops_active)) +
  theme_bw() + 
  #geom_point(aes(size = onset), alpha=.9) + 
  geom_point() +
  geom_smooth(lty = 2, se = FALSE, colour = "firebrick") +
  #theme(plot.title = element_text(face="bold")) +
  scale_y_continuous(limits = c(0, 40)) +
  xlab("Year of Operation") +
  ylab("Active Operations Per Year") +
  ggtitle("The O'Rourke Dataset")

pooled_active <- ggplot(covert_frequency, aes(x = year, y = pooled_ops_active)) +
  theme_bw() + 
  #geom_point(aes(size = onset), alpha=.9) + 
  geom_point() +
  geom_smooth(lty=2, se=FALSE, colour="darkblue") +
  #theme(plot.title = element_text(face="bold")) +
  scale_y_continuous(limits = c(0, 40)) +
  xlab("Year of Operation") +
  #ylab("Active Operations Per Year") +
  ylab("") +
  ggtitle("The Pooled Dataset")

# Plot Grid Guide
# https://cran.r-project.org/web/packages/cowplot/vignettes/plot_grid.html

active_plot <- plot_grid(orke_active, pooled_active)


# save_plot("~/Desktop/active_ops_plot.pdf", active_plot, base_width = 12)

#----- Operations Initiated Per-Year for the O'Rourke and Pooled Samples 

orke_onset <- ggplot(covert_frequency, aes(x = year, y = orourke_ops_initiated)) +
  theme_bw() + 
  #geom_point(aes(size = onset), alpha=.9) + 
  geom_point() +
  geom_smooth(lty = 2, se = FALSE) +
  scale_y_continuous(limits = c(0, 17)) +
  theme(plot.title = element_text(face="bold")) +
  xlab("Year of Operation") +
  ylab("Initiated Per Year") +
  ggtitle("The O'Rourke Dataset")


pooled_onset <- ggplot(covert_frequency, aes(x = year, y = pooled_ops_initiated)) +
  theme_bw() + 
  #geom_point(aes(size = onset), alpha=.9) + 
  geom_point() +
  geom_smooth(lty = 2, se = FALSE) +
  scale_y_continuous(limits = c(0, 17)) +
  theme(plot.title = element_text(face="bold")) +
  xlab("Year of Operation") +
  # ylab("Initiated Per Year") +
  ylab("") +
  ggtitle("The Pooled Dataset")

onset_plot <- plot_grid(orke_onset, pooled_onset)

# save_plot("initiated_per_year_plot.pdf", onset_plot, base_width = 12)


#----- GML MIDs Initiated Per-Year (Includeds Major and Minor MIDs)

minor_onset <- ggplot(covert_frequency, aes(x = year, y = minor_gml_initiated)) +
  theme_bw() + 
  #geom_point(aes(size = onset), alpha=.9) + 
  geom_point() +
  geom_smooth(lty = 2, se = FALSE) +
  scale_y_continuous(limits = c(0, 10), breaks = c(0, 2, 4, 6, 8, 10)) +
  theme() + # plot.title = element_text(face="bold")
  xlab("Year of Operation") +
  ylab("Initiated Per Year") +
  ggtitle("Minor US MIDs")


major_onset <- ggplot(covert_frequency, aes(x = year, y = major_gml_initiated)) +
  theme_bw() + 
  #geom_point(aes(size = onset), alpha=.9) + 
  geom_point() +
  geom_smooth(lty = 2, se = FALSE) +
  scale_y_continuous(limits = c(0, 10), breaks = c(0, 2, 4, 6, 8, 10)) +
  theme() + # plot.title = element_text(face="bold")
  xlab("Year of Operation") +
  # ylab("Initiated Per Year") +
  ylab("") +
  ggtitle("Major US MIDs")

gml_onset_plot <- plot_grid(minor_onset, major_onset)

#--- Saves the plots 

#save_plot("mids_initiated_per_year_plot.pdf", gml_onset_plot, base_width = 12)


#---- Collapses the Data to Make a Count of Active Years per Country


years_active_by_cntry <- dta %>%
  group_by(ccode_current) %>% 
  dplyr::summarise(years_active = sum(orourke_covert_intervention, na.rm = TRUE)) %>% 
  dplyr::rename(ccode = ccode_current)

cshape_1989 <- cshp(date = as.Date("1989-12-31"))

cshape_1989@data = data.frame(cshape_1989@data, 
                              years_active_by_cntry[match(cshape_1989$COWCODE,
                                                          years_active_by_cntry$ccode),])

cshape_1989$operations[is.na(cshape_1989$years_active)] <- 0 

#png('~/Desktop/covert_map.png', units = "in", width = 16, height = 8, res = 500) 


shades = shading(breaks = c(1, 5, 10, 15, 19), cols = brewer.pal(6, "Greys"))
choropleth(cshape_1989, cshape_1989$years_active, shading = shades, border="white", lwd = .25) 
legend(-175, 5, c("0", "1 - 5", "5 - 10", "10 - 15", "15 - 19", "20 - 22"),
       pch = 19, col = paste(brewer.pal(6, "Greys"), sep = ""), 
       cex = 1.6, box.col = 0, title = "Years Active")

# dev.off()


#-----  Plots of Yearly Public Opinion By Quarter 

ggplot(data = us_domestic_data, aes(x = year)) + 
  geom_line(aes(y = approval_q1, colour = "Quarter 1")) + 
  geom_line(aes(y = approval_q2, colour = "Quarter 2")) + 
  geom_line(aes(y = approval_q3, colour = "Quarter 3")) + 
  geom_line(aes(y = approval_q4, colour = "Quarter 4")) + 
  geom_line(aes(y = rand_approv, colour = "Random")) + 
  ylim(c(0, 100)) +
  theme_bw() +
  ylab(label = "Quarterly Public Approval") +  xlab("Year") +
  labs(colour = "Period") 

# ggtitle("Comparison of Quarterly Public Opinion by Year") 

# ggsave("approval_plot.pdf", width = 27, height = 16, units = "cm")


#----- Plots the Bivariate Probit Marginal Effects and Creates a Latex Table
# Note: The BVP Models are run in STATA 

# Prevents R from using scientific notation when displaying the results 
options(scipen = 999)

#----- BVP Onset ME Visualizations 
#----- Creates Subsets of the Data 

dta_me_onset <- bvp_onset_dta %>% filter(Value == "Marginal Effect")
dta_ci_onset <- bvp_onset_dta %>% filter(Value != "Marginal Effect")


#---- Converts all the ME Vectors to Numeric 

dta_me_onset <- dta_me_onset %>% 
  mutate_at(vars(c(3:6)), as.numeric) %>% 
  mutate_if(is.numeric, ~round(., 5))  

glimpse(dta_me_onset)


#----- Merges the Dataframes back together 

dta_clean_onset <- rbind(dta_me_onset, dta_ci_onset) %>% arrange(`Model  (Y1 & Y2)`)

#---- Creates a Latex Table from the Tidy AME Data frame 

kable(dta_clean_onset, "latex", booktabs = T, align=rep('c', 6),
      caption = "The Marginal Effect of Divided Government on the Joint Probability
      of Initiating Covert Operations and Military Force") %>%
  kable_styling(latex_options = c("striped", "hold_position", "scale_down"))


#------ Cleans up the ME Results so that they can be Plotted  

onset_me_tidy <- dta_clean_onset %>% 
  dplyr::rename(measure = Value) %>% 
  gather(key = "key", value = "value", c(3:6)) %>% 
  clean_names() %>% 
  filter(measure == "Marginal Effect") %>% 
  mutate(value = as.numeric(value)) %>% 
  dplyr::select(-measure)

glimpse(onset_me_tidy)

#------ Cleans the Standard errors of the MEs 
onset_se_tidy <- dta_clean_onset %>% 
  dplyr::rename(measure = Value) %>% 
  gather(key = "key", value = "value", c(3:6)) %>% 
  clean_names() %>% 
  filter(measure != "Marginal Effect") %>% 
  mutate(value = gsub("[()]", "", value)) %>%  # removes the parentheses 
  separate(value, c("lower", "upper"), ",") %>% 
  mutate(lower = as.numeric(lower),
         upper = as.numeric(upper)) %>% 
  dplyr::select(-measure)

glimpse(onset_se_tidy)


#----- Combines the Datasets into a Single Dataframe

onset_final_tidy <- left_join(onset_me_tidy, onset_se_tidy)

#---- Creates a Coefficeint Plot from the Tidy Joint Marginal Probability Data

onset_me_plot <- ggplot(onset_final_tidy,
                        aes(x = value, y = key, xmin = lower, xmax = upper,
                            height = 0)) + 
  geom_point() +
  geom_errorbarh() +
  geom_vline(xintercept = 0, color = "gray") +
  facet_wrap(~model_y1_y2) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = "none") +
  ylab("") + xlab("The Marginal Effect of Divided Government") #+
#labs(title = "The Marginal Effects of Divided Government",
#     subtitle = "On the Joint Likelihood of Covert Operations and Military Force")


#setwd("~/Desktop")
#cowplot::save_plot("onset_me_plot.pdf", onset_me_plot, base_width = 7.5, base_height = 3.5)

#----- BVP Ongoing ME Visualizations 
#----- Creates Subsets of the Data 

dta_me_ongoing <- bvp_ongoing_dta %>% filter(Value == "Marginal Effect")
dta_ci_ongoing <- bvp_ongoing_dta %>% filter(Value != "Marginal Effect")


#---- Converts all the ME Vectors to Numeric 

dta_me_ongoing <- dta_me_ongoing %>% 
  mutate_at(vars(c(3:6)), as.numeric) %>% 
  mutate_if(is.numeric, ~round(., 5))  

glimpse(dta_me_ongoing)


#----- Merges the Dataframes back together 

dta_clean_ongoing <- rbind(dta_me_ongoing, dta_ci_ongoing) %>% arrange(`Model  (Y1 & Y2)`)

#---- Creates a Latex Table from the Tidy AME Data frame 

kable(dta_clean_ongoing, "latex", booktabs = T, align=rep('c', 6),
      caption = "The Marginal Effect of Divided Government on the Joint Probability
      of Covert Operations and Military Force") %>%
  kable_styling(latex_options = c("striped", "hold_position", "scale_down"))



#------ Cleans up the ME Results so that they can be Plotted  

ongoing_me_tidy <- dta_clean_ongoing %>% 
  dplyr::rename(measure = Value) %>% 
  gather(key = "key", value = "value", c(3:6)) %>% 
  clean_names() %>% 
  filter(measure == "Marginal Effect") %>% 
  mutate(value = as.numeric(value)) %>% 
  dplyr::select(-measure)

glimpse(ongoing_me_tidy)

#------ Cleans the Standard errors of the MEs 
ongoing_se_tidy <- dta_clean_ongoing %>% 
  dplyr::rename(measure = Value) %>% 
  gather(key = "key", value = "value", c(3:6)) %>% 
  clean_names() %>% 
  filter(measure != "Marginal Effect") %>% 
  mutate(value = gsub("[()]", "", value)) %>%  # removes the parentheses 
  separate(value, c("lower", "upper"), ",") %>% 
  mutate(lower = as.numeric(lower),
         upper = as.numeric(upper)) %>% 
  dplyr::select(-measure)

glimpse(ongoing_se_tidy)


#----- Combines the Datasets into a Single Dataframe

ongoing_final_tidy <- left_join(ongoing_me_tidy, ongoing_se_tidy)

#---- Creates a Coefficeint Plot from the Tidy Joint Marginal Probability Data

ongoing_me_plot <- ggplot(ongoing_final_tidy,
                          aes(x = value, y = key, xmin = lower, xmax = upper,
                              height = 0)) + 
  geom_point() +
  geom_errorbarh() +
  geom_vline(xintercept = 0, color = "gray") +
  facet_wrap(~model_y1_y2) +
  theme_bw(base_family = "Helvetica") +
  theme(legend.position = "none") +
  ylab("") + xlab("The Marginal Effect of Divided Government") #+
# labs(title = "The Marginal Effects of Divided Government",
#     subtitle = "On the Joint Likelihood of Covert Operations and Military Force")

# cowplot::save_plot("ongoing_me_plot.pdf", ongoing_me_plot, base_width = 7.5, base_height = 3.5)






